########################################################################
########################################################################
# Replication code for:
#
# @article{CansunarFountains,
# title={Distributional Consequences of Philanthropic Contributions to Public Goods: Self-Serving Elite in Ottoman Istanbul},
# author={Cansunar, Asli},
# journal={the Journal of Politics},
# year={2021},
# }
#
# Comments, questions: asli.cansunar@nuffield.ox.ac.uk
########################################################################
########################################################################

rm(list=ls())

#### 0. Package Dependencies ####
if(!require(pacman)){install.packages("pacman");require(pacman)}
#Install/load tons of packages
p_load(plyr,stargazer, MASS, CBPS, ggplot2, sandwich, 
       rgdal, proj4, maptools, spdep, visreg, lmtest, 
       xtable, CBPS, AER, ggpredict,readstata13,ggmap,splitstackshape,
       dplyr, readr, githubinstall,remotes,sjPlot,sjmisc,lmtest,
       sandwich,multiwayvcov,mgcv,maps,dummies,sf,rspatial,dismo,
       sp,mapview,osmdata,tmap,lwgeom,deldir,voronoi,glmmfields,
       texreg,data.table,coefplot,mfx,spatialreg, pscl,gridExtra,
       huxtable, broom, officer, flextable,readxl)


## set working directory
setwd("###SET YOUR WORKING DIRECTORY")

#### 1. Data ####
newdata <- read.dta13("fountains.dta")

complete.cases(newdata)

#### 2. Spatial Weights ####
df<-subset(newdata, select=c("latitude", "longitude"))
voronoi<-voronoi(df)
plot(voronoi)
proj4string(voronoi)<- CRS("+proj=longlat +datum=WGS84")
s_poly <- SpatialPolygonsDataFrame(voronoi, newdata)
W_nb <- poly2nb(voronoi, queen = T)

## Visualize queen neighbors
plot(voronoi, col="gray", border="blue")
plot(W_nb, df, add=TRUE, pch=20)

## Compute spatial weights
W_cont_matW <- nb2listw(W_nb, style = "W", zero.policy = T)
summary(W_nb)
plot(W_nb,coordinates(df))

#### Figure 1 ####

postscript("Figure1.eps",   width=8, height=5)
p1<-ggplot(df, aes(x=newdata$number)) + geom_histogram(color="black", fill="gray",binwidth=1)+
      labs(x="Number of fountains", y = "Frequency")+  theme_classic()
p2<-ggplot(df, aes(x=newdata$water_per_day)) + geom_histogram(color="black", fill="gray",binwidth=100)+
  labs(x="Cubic meters of water per day", y = "Frequency")+  theme_classic()
grid.arrange(p1, p2, nrow = 1)
dev.off()

#### Figure 2 ####
register_google(key = "###ENTER YOUR OWN API###")
newdata$lon<-newdata$longitude
newdata$lat<-newdata$latitude
mean.longitude <- mean(newdata$longitude)
mean.latitude <- mean(newdata$latitude)
n.map <- get_map(location = c(mean.longitude, mean.latitude), zoom = 13, scale = 2,source = "stamen", maptype = "toner")
dfexp <- expandRows(newdata, "number")
data.muslim <- subset(newdata, NonMuslim==0)
data.nmuslim <- subset(newdata, NonMuslim==1)
data.elite <- subset(newdata, elite_house==1)
data.commoner <- subset(newdata, elite_house==0)

ggsave("Figure2.eps", device=cairo_ps)
ggmap(n.map) +  
  aes(x=lon, y=lat) + 
  stat_density_2d(aes(fill= stat(nlevel)), dfexp, geom="polygon",alpha = .6, bins=20) + 
  scale_fill_gradientn(colors = brewer.pal(7, "YlOrRd"))+
  geom_point(data=data.elite, position="jitter", alpha=0.8, size=1, shape=15) +
  theme(
    legend.position='none'
  )
dev.off()


#### Figure 3 ####
ggsave("Figure3.eps", device=cairo_ps)
ggmap(n.map) +  
  aes(x=lon, y=lat)+ 
  stat_density_2d(aes(fill= stat(nlevel)), dfexp, geom="polygon",alpha = .6, bins=20) + 
  scale_fill_gradientn(colors = brewer.pal(7, "YlOrRd"))+
  geom_point(data=data.nmuslim, position="jitter", alpha=0.8, size=0.8, shape=17) +
  theme(
    legend.position='none'
  )
dev.off()


#### Table 1 ####

m1<-glm.nb(number~ NonMuslim + No_household, data=newdata)
mirr1<-negbinirr(number~ NonMuslim + No_household, data=newdata)
summary(m1)
resm1 <- m1$residuals
moran.mc(resm1, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.05 indicates no spatial autocorrelation and no need for spatial filtering

m2<-glm.nb(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
mirr2<-negbinirr(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
summary(m2)
resm2 <- m2$residuals
moran.mc(resm2, listw=W_cont_matW, zero.policy=T,nsim=10000)

#Moran I = 0.04 indicates no spatial autocorrelation and no need for spatial filtering

m3<-glm.nb(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
mirr3<-negbinirr(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
summary(m3)
resm3 <- m3$residuals
moran.mc(resm3, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.04 indicates no spatial autocorrelation and no need for spatial filtering

m4<-glm.nb(formula=number~ elite_house+ No_household, data=newdata)
mirr4<-negbinirr(formula=number~ elite_house+ No_household, data=newdata)
summary(m4)
resm4 <- m4$residuals
moran.mc(resm4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I =0.05 indicates no spatial autocorrelation and no need for spatial filtering

m5<-glm.nb(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
mirr5<-negbinirr(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
summary(m5)
resm5 <- m5$residuals
moran.mc(resm5, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.05 indicates no spatial autocorrelation and no need for spatial filtering

m6<-glm.nb(number~ elite_house+ No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
mirr6<-negbinirr(number~ elite_house+ No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
summary(m6)
resm6 <- m6$residuals
moran.mc(resm6, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I =0.05 indicates no spatial autocorrelation and no need for spatial filtering

m7<-glm.nb(number~ as.factor(NonMuslim) + as.factor(elite_house) + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata, link=log)
summary(m7)
mirr7<-negbinirr(number~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata)
resm7 <- m7$residuals
moran.mc(resm7, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.02 indicates no spatial autocorrelation and no need for spatial filtering

#Table 1 in the main text
texreg(list(mirr1,mirr2,mirr3,mirr4,mirr5,mirr6,mirr7), stars=c(0.01, 0.05,0.1), digits=3,reorder.coef = c(1, 3, 2), file="table1.tex",omit.coef="(Elevation)|(WalkingDistanceKilometers)|(StraightLineDistanceKilomete)|(Eyup)",custom.coef.names = c("Non-Muslim Neighborhood", "No. of households","Elite neighborhood"))

#Full Table - Table A2 in the appendix
texreg(list(mirr1,mirr2,mirr3,mirr4,mirr5,mirr6,mirr7), stars=c(0.01, 0.05,0.1), digits=3,reorder.coef = c(1, 7, 2,3,4,5,6), file="tableA2.tex",custom.coef.names = c("Non-Muslim Neighborhood", "No. of households","Elevation","Distance to the palace","Distance fromt the coast", "Eyup","Elite neighborhood"))

#### Table 2 ####

ols1<-lm(water_per_day~ NonMuslim + No_household, data=newdata)
summary(ols1)
resols1 <- ols1$residuals
moran.test(resols1, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.17 indicates mild spatial autocorrelation
Eigenols1<-SpatialFiltering(water_per_day ~ NonMuslim + No_household, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
ols1Eig<-lm(water_per_day~ NonMuslim + No_household+fitted(Eigenols1), data=newdata)
summary(ols1Eig)


ols2<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
summary(ols2)
resols2 <- ols2$residuals
moran.test(resols2, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15  indicates mild spatial autocorrelation
Eigenols2<-SpatialFiltering(water_per_day ~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
ols2Eig<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(Eigenols2), data=newdata)
summary(ols2Eig)

ols3<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
summary(ols3)
resols3 <- ols3$residuals
moran.test(resols3, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 indicates mild spatial autocorrelation
Eigenols3<-SpatialFiltering(water_per_day ~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
ols3Eig<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup+fitted(Eigenols3), data=newdata)
summary(ols3Eig)


ols4<-lm(water_per_day~ elite_house + No_household, data=newdata)
summary(ols4)
resols4 <- ols4$residuals
moran.test(resols4, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.2 indicates mild spatial autocorrelation
Eigenols4<-SpatialFiltering(water_per_day ~ elite_house + No_household, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
ols4Eig<-lm(water_per_day~ elite_house + No_household+fitted(Eigenols4), data=newdata)
summary(ols4Eig)


ols5<-lm(water_per_day~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
summary(ols5)
resols5 <- ols5$residuals
moran.test(resols5, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.18 indicates mild spatial autocorrelation
Eigenols5<-SpatialFiltering(water_per_day ~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
ols5Eig<-lm(water_per_day~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(Eigenols5), data=newdata)
summary(ols5Eig)

ols6<-lm(water_per_day~ elite_house+ No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
summary(ols6)
resols6 <- ols6$residuals
moran.test(resols6, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.17  indicates mild spatial autocorrelation
Eigenols6<-SpatialFiltering(water_per_day ~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
ols6Eig<-lm(water_per_day~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup+fitted(Eigenols6), data=newdata)
summary(ols6Eig)
moran.plot(resols6, 
           W_cont_matW)
spplot(voronoi,zcol="number$newdata",col="transparent")


ols7<-lm(water_per_day~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata)
summary(ols7)
resols7 <- ols7$residuals
moran.test(resols7, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 with p=7.161e-07 indicates mild spatial autocorrelation
Eigenols7<-SpatialFiltering(water_per_day ~ NonMuslim+elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
ols7Eig<-lm(water_per_day~ NonMuslim+ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup+fitted(Eigenols7), data=newdata)
summary(ols7Eig)

stargazer(ols1Eig,ols2Eig,ols3Eig,ols4Eig,ols5Eig,ols6Eig,ols7Eig,digits=2,omit=c("Elevation","WalkingDistanceKilometers","StraightLineDistanceKilomete","Eyup","fitted"), covariate.labels = c("Non-Muslim neighborhood","Elite neighborhood", "No. of households","Constant"), omit.stat=c("rsq","ser","f"), add.lines=list(c("Geographic Controls", "No", "No", "Yes","No", "No", "Yes", "Yes"), c("District FE", "No", "No", "Yes", "No", "No", "Yes", "Yes"), c("Moran Eigenvectors", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")), column.labels=c("Model 1", "Model 2", "Model 3", "Model 4","Model 5","Model 6","Model 7"), star.cutoffs=c(0.1, 0.05, 0.01), star.char=c("*", "**", "***"), notes = c("* p<0.1; ** p<0.05; *** p<0.01"), notes.append = F,out="table2.tex")
stargazer(ols1Eig,ols2Eig,ols3Eig,ols4Eig,ols5Eig,ols6Eig,ols7Eig,digits=2, omit=c("fitted"),omit.stat=c("rsq","ser","f"), column.labels=c("Model 1", "Model 2", "Model 3", "Model 4","Model 5","Model 6","Model 7"),covariate.labels = c("Non-Muslim Neighborhood", "Elite neighborhood","No. of households","Elevation","Distance to the palace","Distance fromt the coast", "Eyup","Constant"), star.cutoffs=c(0.1, 0.05, 0.01), star.char=c("*", "**", "***"), notes = c("* p<0.1; ** p<0.05; *** p<0.01"),out="tableA3.tex")

#### Table 3 ####

ch1<-hurdle(number~ NonMuslim  + No_household, data=newdata,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
summary(ch1)
resch1 <- ch1$residuals
moran.mc(resch1, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.05 indicates no spatial autocorrelation

ch2<-hurdle(number~ NonMuslim  + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch2<-residuals(ch2)
moran.mc(resch2, listw=W_cont_matW, zero.policy=T,nsim=10000)
summary(ch2)
#Moran I statistic of 0.04 indicates no spatial autocorrelation

ch3<-hurdle(number~ NonMuslim  + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch3<-residuals(ch3)
moran.mc(resch3, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.03 indicates no spatial autocorrelation


ch4<-hurdle(number~ elite_house + No_household, data=newdata,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
summary(ch4)
resch4 <- ch4$residuals
moran.mc(resch4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.05 indicates no spatial autocorrelation


ch5<-hurdle(number~ elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch5<-residuals(ch2)
moran.mc(resch5, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.04 indicates no spatial autocorrelation
summary(ch5)

ch6<-hurdle(number~ elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch6<-residuals(ch6)
moran.mc(resch6, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.04 indicates no spatial autocorrelation
summary(ch6)

ch7<-hurdle(number~ NonMuslim+elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch7<-residuals(ch7)
moran.mc(resch7, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.02 indicates no spatial autocorrelation
summary(ch7)


texreg(list(ch1,ch2,ch3,ch4,ch5,ch6,ch7), stars=c(0.01, 0.05,0.1), digits=3,booktabs=TRUE,file="table3.tex", omit.coef=("(WalkingDistanceKilometers)|(Elevation)|(Eyup)|(StraightLineDistanceKilomete)|(Intercept)|(Log(theta))"))


###Figure 4 and 5
cp <- glm(factor(number > 0) ~ NonMuslim + elite_house+No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup,
          data = newdata, family = binomial(link = "probit"))
v<-visreg(cp, "No_household",by="elite_house",scale = "response", xlab="Number of households", ylab="Probability of a fountain",overlay=TRUE,legend=FALSE,line=list(lty=1:2, col="black"), fill=list(col=grey(c(0.2,0.8), alpha=0.4)))

pdf("Figure5.pdf",   width=6, height=5)
par(mar = c(8,4,4,4))
plot(v, overlay=TRUE, legend=FALSE,
     ylab="Predicted probability of a fountain",
     xlab="Number of households",
     font.lab=1,
     yaxs="i",
     ylim=c(0.5,1),
     xlim=c(0,300),
     cex.lab=1,
     cex.lab=1,line=list(lty=1:2, col="black"),fill=list(col=grey(c(0.2,0.8), alpha=0.4)))
legend("bottom", inset=c(0,-0.6), ,lty=1:2,legend=c("Commoner","Elite"),xpd=TRUE,horiz = TRUE)
dev.off()


v2<-visreg(cp, "No_household",by="NonMuslim",scale = "response", xlab="Number of households", 
           ylab="Probability of a fountain",overlay=TRUE,legend=FALSE,
           line=list(lty=1:2, col="black"), fill=list(col=grey(c(0.2,0.8), alpha=0.4)))


pdf("Figure4.pdf",   width=6, height=5)
par(mar = c(8,4,4,4))
plot(v2, overlay=TRUE, legend=FALSE,
     ylab="Predicted probability of a fountain",
     xlab="Number of households",
     font.lab=1,
     yaxs="i",
     ylim=c(0.3,1),
     xlim=c(0,300),
     cex.lab=1,line=list(lty=1:2, col="black"),fill=list(col=grey(c(0.2,0.8), alpha=0.4)))
legend("bottom", inset=c(0,-0.6), lty=1:2,legend=c("Muslim","Non-Muslim"),xpd=TRUE,horiz = TRUE)
dev.off()


#### Table 4 ####

olsLC1<-lm(least_cost~ NonMuslim + No_household, data=newdata)
summary(olsLC)
resolsLC1 <- olsLC1$residuals
moran.test(resolsLC1, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.18 with p-value = 1.034e-08 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC1<-SpatialFiltering(least_cost~ NonMuslim+ No_household, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
olsLCEig1<-lm(least_cost~ NonMuslim + No_household +fitted(EigenolsLC1), data=newdata)
summary(olsLCEig)

olsLC2<-lm(least_cost~ NonMuslim  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata)
summary(olsLC2)
resolsLC2 <- olsLC2$residuals
moran.test(resolsLC2, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.16 with p-value = 2.395e-07 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC2<-SpatialFiltering(least_cost~ NonMuslim  + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
olsLCEig2<-lm(least_cost~ NonMuslim  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(EigenolsLC2), data=newdata)
summary(olsLCEig)

olsLC3<-lm(least_cost~ NonMuslim  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata)
summary(olsLC3)
resolsLC3 <- olsLC3$residuals
moran.test(resolsLC3, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.16 with p-value = 2.395e-07 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC3<-SpatialFiltering(least_cost~ NonMuslim + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
olsLCEig3<-lm(least_cost~ NonMuslim  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup+fitted(EigenolsLC3), data=newdata)


olsLC4<-lm(least_cost~ elite_house + No_household, data=newdata)
resolsLC4 <- olsLC4$residuals
moran.test(resolsLC4, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.17 with p-value = 6.381e-08 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC4<-SpatialFiltering(least_cost~ elite_house+ No_household, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
olsLCEig4<-lm(least_cost~ elite_house + No_household +fitted(EigenolsLC4), data=newdata)


olsLC5<-lm(least_cost~ elite_house  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata)
summary(olsLC5)
resolsLC5 <- olsLC5$residuals
moran.test(resolsLC5, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 with p-value = 1.404e-06 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC5<-SpatialFiltering(least_cost~ elite_house  + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
olsLCEig5<-lm(least_cost~ elite_house  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(EigenolsLC5), data=newdata)


olsLC6<-lm(least_cost~ elite_house  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata)
summary(olsLC6)
resolsLC6 <- olsLC6$residuals
moran.test(resolsLC6, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 with p-value = 1.404e-06 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC6<-SpatialFiltering(least_cost~ elite_house + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
olsLCEig6<-lm(least_cost~ elite_house  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup+fitted(EigenolsLC6), data=newdata)

olsLC7<-lm(least_cost~ NonMuslim + elite_house + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata)
summary(olsLC7)
resolsLC7 <- olsLC7$residuals
moran.test(resolsLC7, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 with p-value = 1.412e-06 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC7<-SpatialFiltering(least_cost~ NonMuslim + elite_house + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=newdata, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #
olsLCEig7<-lm(least_cost~ NonMuslim + elite_house + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup+fitted(EigenolsLC7), data=newdata)

stargazer(olsLCEig1,olsLCEig2,olsLCEig3,olsLCEig4,olsLCEig5,olsLCEig6,olsLCEig7,digits=2,omit=c("WalkingDistanceKilometers","StraightLineDistanceKilomete","Eyup","fitted"), covariate.labels = c("Non-Muslim neighborhood","Elite neighborhood", "No. of households","Constant"), omit.stat=c("rsq","ser","f"), add.lines=list(c("Geographic Controls", "No", "Yes", "Yes","No", "Yes", "Yes", "Yes"), c("District FE", "No", "No", "Yes", "No", "No", "Yes", "Yes"), c("Moran Eigenvectors", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")), column.labels=c("Model 1", "Model 2", "Model 3", "Model 4","Model 5","Model 6","Model 7"), star.cutoffs=c(0.1, 0.05, 0.01), star.char=c("*", "**", "***"), notes = c("* p<0.1; ** p<0.05; *** p<0.01"), notes.append = F,out="Table4.tex")


######################################################################
#### 4. Appendix ####

#### Figure A1 ####
data.nowater <- subset(newdata, Water_existence==0)
data.water<- subset(newdata, Water_existence==1)
map <- get_googlemap(center = c(mean.longitude, mean.latitude), maptype = "terrain", source = "google", zoom = 13, color = "bw")
basemap_attributes <- attributes(map)
basemap_transparent <- matrix(adjustcolor(map, 
                                          alpha.f = 0.5), 
                              nrow = nrow(map))
attributes(basemap_transparent) <- basemap_attributes
longitude<-28.983551
latitude<-41.011625
palace <- data.frame(longitude,latitude)

pdf("FigureA1.pdf",   width=7, height=5)
ggmap(basemap_transparent) + geom_point(
  aes(x = longitude, y = latitude),
  data = newdata, colour = "red", position="jitter", alpha=0.8,shape=19, size=0.9)+
  geom_point(
    aes(x = longitude, y = latitude),
    data = palace, colour = "purple", alpha=0.8,shape=18, size=5)
dev.off()

#### Figure A2 ####
pdf("FigureA2.pdf",   width=7, height=5)
ggmap(basemap_transparent) + geom_point(
  aes(x = longitude, y = latitude),
  data = data.nmuslim, colour = "red", position="jitter", alpha=0.8, size=1.2,shape=17)+
  geom_point(
    aes(x = longitude, y = latitude),
    data = palace, colour = "purple", alpha=0.8,shape=18, size=5)
dev.off()


#### Figure A3 ####
church<- read_excel("nmuslim.xlsx")

pdf("FigureA3.pdf",   width=7, height=5)
ggmap(basemap_transparent) + geom_point(
  aes(x = longitude, y = latitude),
  data = data.nmuslim, colour = "red", position="jitter", alpha=0.8, size=1,shape=17)+geom_point(
    aes(x = Lon, y = Lat),
    data = church, colour = "blue", position="jitter", alpha=0.5, size=1)
dev.off()

#### Figure A4 ####
pdf("FigureA4.pdf",   width=7, height=5)
ggmap(basemap_transparent) + geom_point(
  aes(x = longitude, y = latitude),
  data = data.elite, colour = "red", position="jitter", alpha=0.8, size=1.2,shape=19)+
  geom_point(
    aes(x = longitude, y = latitude),
    data = palace, colour = "purple", alpha=0.8,shape=18, size=5)
dev.off()
#### Table A1 ####
stargazer(newdata,mydatastyle="apsr", type = "text",out="TableA1.tex")

#### Table A4 ####


m1<-glm.nb(number~ NonMuslim*No_household, data=newdata)
mirr1<-negbinirr(number~ NonMuslim * No_household, data=newdata)
summary(m1)
resm1 <- m1$residuals
moran.mc(resm1, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.05 indicates no spatial autocorrelation and no need for spatial filtering

m2<-glm.nb(number~ NonMuslim*No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
mirr2<-negbinirr(number~ NonMuslim * No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
summary(m2)
resm2 <- m2$residuals
moran.mc(resm2, listw=W_cont_matW, zero.policy=T,nsim=10000)

#Moran I = 0.04 indicates no spatial autocorrelation and no need for spatial filtering

m3<-glm.nb(number~ NonMuslim* No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
mirr3<-negbinirr(number~ NonMuslim * No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
summary(m3)
resm3 <- m3$residuals
moran.mc(resm3, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.04 indicates no spatial autocorrelation and no need for spatial filtering

m4<-glm.nb(formula=number~ elite_house*No_household, data=newdata)
mirr4<-negbinirr(formula=number~ elite_house* No_household, data=newdata)
summary(m4)
resm4 <- m4$residuals
moran.mc(resm4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I =0.05 indicates no spatial autocorrelation and no need for spatial filtering

m5<-glm.nb(number~ elite_house*No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
mirr5<-negbinirr(number~ elite_house * No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=newdata)
summary(m5)
resm5 <- m5$residuals
moran.mc(resm5, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.05 indicates no spatial autocorrelation and no need for spatial filtering

m6<-glm.nb(number~ elite_house*No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
mirr6<-negbinirr(number~ elite_house*No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=newdata)
summary(m6)
resm6 <- m6$residuals
moran.mc(resm6, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I =0.05 indicates no spatial autocorrelation and no need for spatial filtering


texreg(list(mirr1,mirr2,mirr3,mirr4,mirr5,mirr6), stars=c(0.01, 0.05,0.1), digits=3, file="tableA5.tex")

###Figure A5 ###


v1<-visreg(m3, "No_household",by="NonMuslim", xlab="Number of households", 
           ylab="Predicted number of fountains",overlay=TRUE,legend=FALSE,
           line=list(lty=1:2, col="black"), fill=list(col=grey(c(0.2,0.8), alpha=0.4)), partial=FALSE)

pdf("FigureA5.pdf",   width=6, height=5)
par(mar = c(8,4,4,4))
plot(v1, overlay=TRUE, legend=FALSE,partial=FALSE,
     ylab="Predicted number of fountains",
     xlab="Number of households",
     font.lab=1,
     yaxs="i",
     ylim=c(0,2),
     xlim=c(0,300),
     cex.lab=1,line=list(lty=1:2, col="black"),fill=list(col=grey(c(0.2,0.8), alpha=0.4)))
legend("bottom", inset=c(0,-0.6), lty=1:2,legend=c("Muslim","Non-Muslim"),xpd=TRUE,horiz = TRUE)
dev.off()

v2<-visreg(m6, "No_household",by="elite_house", xlab="Number of households", 
           ylab="Probability of a fountain",overlay=TRUE,legend=FALSE,
           line=list(lty=1:2, col="black"), fill=list(col=grey(c(0.2,0.8), alpha=0.4)),partial=FALSE)

###Figure A6 ###
pdf("FigureA6.pdf",   width=6, height=5)
par(mar = c(8,4,4,4))
plot(v2, overlay=TRUE, legend=FALSE,partial=FALSE,
     ylab="Predicted number of fountains",
     xlab="Number of households",
     font.lab=1,
     yaxs="i",
     ylim=c(0,2),
     xlim=c(0,300),
     cex.lab=1,
     cex.lab=1,line=list(lty=1:2, col="black"),fill=list(col=grey(c(0.2,0.8), alpha=0.4)))
legend("bottom", inset=c(0,-0.6), ,lty=1:2,legend=c("Commoner","Elite"),xpd=TRUE,horiz = TRUE)
dev.off()



### Testing spatial auto-correlation
moran.mc(newdata$elite_house, listw=W_cont_matW, zero.policy=T,nsim=10000)
moran.mc(newdata$NonMuslim, listw=W_cont_matW, zero.policy=T,nsim=10000)
moran.mc(newdata$number, listw=W_cont_matW, zero.policy=T,nsim=10000)
moran.mc(newdata$water_per_day, listw=W_cont_matW, zero.policy=T,nsim=10000)
moran.mc(newdata$least_cost, listw=W_cont_matW, zero.policy=T,nsim=10000)


### Figure A8
pdf("figure_fountainmoran.pdf",   width=7, height=5)
moran.plot(newdata$number, listw=W_cont_matW, xlab="Number of fountains", ylab="Spatially lagged number of fountains")
dev.off()

### Figure A9
pdf("figure_watermoran.pdf",   width=7, height=5)
moran.plot(newdata$water_per_day, listw=W_cont_matW, xlab="Water per day", ylab="Spatially lagged water per day")
dev.off()

### Figure A10
pdf("figure_flcmoran.pdf",   width=7, height=5)
moran.plot(newdata$least_cost, listw=W_cont_matW, xlab="Least cost", ylab="Spatially lagged least cost")
dev.off()



#### Table A5 ####
intramuros <- subset(newdata, Eyup==0)
df<-subset(intramuros, select=c("latitude", "longitude"))
voronoi<-voronoi(df)
plot(voronoi)
proj4string(voronoi)<- CRS("+proj=longlat +datum=WGS84")

W_nb <- poly2nb(voronoi, queen = T)
##visualize queen neighbors
plot(voronoi, col="gray", border="blue")
plot(W_nb, df, add=TRUE, pch=20)
#Compute spatial weights:  A spatial weights matrix reflects the intensity of the geographic relationship between observations
W_cont_matW <- nb2listw(W_nb, style = "W", zero.policy = T)
summary(W_nb)
plot(W_nb,coordinates(df))

m1<-glm.nb(number~ NonMuslim + No_household, data=intramuros)
mirr1<-negbinirr(number~ NonMuslim + No_household, data=intramuros)
summary(m1)
resm1 <- m1$residuals
moran.mc(resm1, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.03 indicates no spatial autocorrelation and no need for spatial filtering

m2<-glm.nb(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
mirr2<-negbinirr(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(m2)
resm2 <- m2$residuals
moran.mc(resm2, listw=W_cont_matW, zero.policy=T,nsim=10000)

#Moran I = 0.03 indicates no spatial autocorrelation and no need for spatial filtering


m4<-glm.nb(formula=number~ elite_house+ No_household, data=intramuros)
mirr4<-negbinirr(formula=number~ elite_house+ No_household, data=intramuros)
summary(m4)
resm4 <- m4$residuals
moran.mc(resm4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I =0.04 indicates no spatial autocorrelation and no need for spatial filtering

m5<-glm.nb(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
mirr5<-negbinirr(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(m5)
resm5 <- m5$residuals
moran.mc(resm5, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.04 indicates no spatial autocorrelation and no need for spatial filtering


m7<-glm.nb(number~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
mirr7<-negbinirr(number~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=intramuros)
resm7 <- m7$residuals
moran.mc(resm7, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.02 indicates no spatial autocorrelation and no need for spatial filtering

texreg(list(mirr1,mirr2,mirr4,mirr5,mirr7), stars=c(0.01, 0.05,0.1), digits=3,reorder.coef = c(1, 3, 2), file="nbreg_intramuros.tex",omit.coef="(Elevation)|(WalkingDistanceKilometers)|(StraightLineDistanceKilomete)|(Eyup)",custom.coef.names = c("Non-Muslim Neighborhood", "No. of households","Elite neighborhood"),out="TableA5.tex")

#### Table A6 ####
ols1<-lm(water_per_day~ NonMuslim + No_household, data=intramuros)
summary(ols1)
resols1 <- ols1$residuals
moran.test(resols1, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.16 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols1<-SpatialFiltering(water_per_day ~ NonMuslim + No_household, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols1Eig<-lm(water_per_day~ NonMuslim + No_household+fitted(Eigenols1), data=intramuros)
summary(ols1Eig)


ols2<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(ols2)
resols2 <- ols2$residuals
moran.test(resols2, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15  indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols2<-SpatialFiltering(water_per_day ~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols2Eig<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(Eigenols2), data=intramuros)
summary(ols2Eig)



ols4<-lm(water_per_day~ elite_house + No_household, data=intramuros)
summary(ols4)
resols4 <- ols4$residuals
moran.test(resols4, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.2 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols4<-SpatialFiltering(water_per_day ~ elite_house + No_household, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols4Eig<-lm(water_per_day~ elite_house + No_household+fitted(Eigenols4), data=intramuros)
summary(ols4Eig)


ols5<-lm(water_per_day~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(ols5)
resols5 <- ols5$residuals
moran.test(resols5, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.18 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols5<-SpatialFiltering(water_per_day ~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols5Eig<-lm(water_per_day~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(Eigenols5), data=intramuros)
summary(ols5Eig)


ols7<-lm(water_per_day~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(ols7)
resols7 <- ols7$residuals
moran.test(resols7, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 with p=7.161e-07 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols7<-SpatialFiltering(water_per_day ~ NonMuslim+elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols7Eig<-lm(water_per_day~ NonMuslim+ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup+fitted(Eigenols7), data=intramuros)
summary(ols7Eig)

stargazer(ols1Eig,ols2Eig,ols4Eig,ols5Eig,ols7Eig,digits=2,omit=c("Elevation","WalkingDistanceKilometers","StraightLineDistanceKilomete","fitted"), covariate.labels = c("Non-Muslim neighborhood","Elite neighborhood", "No. of households","Constant"), omit.stat=c("rsq","ser","f"), add.lines=list(c("Geographic Controls", "No", "Yes","No", "Yes", "Yes"), c("Moran Eigenvectors", "Yes", "Yes", "Yes", "Yes", "Yes")), column.labels=c("Model 1", "Model 2", "Model 3", "Model 4","Model 5"), star.cutoffs=c(0.1, 0.05, 0.01), star.char=c("*", "**", "***"), notes = c("* p<0.1; ** p<0.05; *** p<0.01"), notes.append = F,out="TableA6.tex")


#### Table A7 ####
ch1<-hurdle(number~ NonMuslim  + No_household, data=intramuros,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
summary(ch1)
resch1 <- ch1$residuals
moran.mc(resch1, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.05 indicates no spatial autocorrelation

ch2<-hurdle(number~ NonMuslim  + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch2<-residuals(ch2)
moran.mc(resch2, listw=W_cont_matW, zero.policy=T,nsim=10000)
summary(ch2)
#Moran I statistic of 0.04 indicates no spatial autocorrelation

ch4<-hurdle(number~ elite_house + No_household, data=intramuros,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
summary(ch4)
resch4 <- ch4$residuals
moran.mc(resch4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.05 indicates no spatial autocorrelation


ch5<-hurdle(number~ elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch5<-residuals(ch2)
moran.mc(resch5, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.04 indicates no spatial autocorrelation
summary(ch5)


ch7<-hurdle(number~ NonMuslim+elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch7<-residuals(ch7)
moran.mc(resch7, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.02 indicates no spatial autocorrelation
summary(ch7)

texreg(list(ch1,ch2,ch4,ch5,ch7), stars=c(0.01, 0.05,0.1), digits=3,file="TableA7.tex")


#### Table A8 ####
olsLC1<-lm(least_cost~ NonMuslim + No_household, data=intramuros)
summary(olsLC)
resolsLC1 <- olsLC1$residuals
moran.test(resolsLC1, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.18 with p-value = 1.034e-08 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC1<-SpatialFiltering(least_cost~ NonMuslim+ No_household, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
olsLCEig1<-lm(least_cost~ NonMuslim + No_household +fitted(EigenolsLC1), data=intramuros)
summary(olsLCEig)

olsLC2<-lm(least_cost~ NonMuslim  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(olsLC2)
resolsLC2 <- olsLC2$residuals
moran.test(resolsLC2, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.16 with p-value = 2.395e-07 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC2<-SpatialFiltering(least_cost~ NonMuslim  + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
olsLCEig2<-lm(least_cost~ NonMuslim  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(EigenolsLC2), data=intramuros)
summary(olsLCEig)




olsLC4<-lm(least_cost~ elite_house + No_household, data=intramuros)
resolsLC4 <- olsLC4$residuals
moran.test(resolsLC4, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.17 with p-value = 6.381e-08 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC4<-SpatialFiltering(least_cost~ elite_house+ No_household, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
olsLCEig4<-lm(least_cost~ elite_house + No_household +fitted(EigenolsLC4), data=intramuros)


olsLC5<-lm(least_cost~ elite_house  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(olsLC5)
resolsLC5 <- olsLC5$residuals
moran.test(resolsLC5, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 with p-value = 1.404e-06 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC5<-SpatialFiltering(least_cost~ elite_house  + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
olsLCEig5<-lm(least_cost~ elite_house  + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(EigenolsLC5), data=intramuros)


olsLC7<-lm(least_cost~ NonMuslim + elite_house + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros)
summary(olsLC7)
resolsLC7 <- olsLC7$residuals
moran.test(resolsLC7, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.15 with p-value = 1.412e-06 indicates mild spatial autocorrelation
old <- Sys.time()
EigenolsLC7<-SpatialFiltering(least_cost~ NonMuslim + elite_house + No_household + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=intramuros, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
olsLCEig7<-lm(least_cost~ NonMuslim + elite_house + No_household  + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(EigenolsLC7), data=intramuros)

stargazer(olsLCEig1,olsLCEig2,olsLCEig4,olsLCEig5,olsLCEig7,digits=2,omit=c("WalkingDistanceKilometers","StraightLineDistanceKilomete","fitted"), covariate.labels = c("Non-Muslim neighborhood","Elite neighborhood", "No. of households","Constant"), omit.stat=c("rsq","ser","f"), add.lines=list(c("Geographic Controls", "No", "Yes","No", "Yes", "Yes"), c("Moran Eigenvectors", "Yes", "Yes", "Yes", "Yes", "Yes")), column.labels=c("Model 1", "Model 2", "Model 3", "Model 4","Model 5"), star.cutoffs=c(0.1, 0.05, 0.01), star.char=c("*", "**", "***"), notes = c("* p<0.1; ** p<0.05; *** p<0.01"), notes.append = F,out="tableA8.tex")

### Table A10 ###


df<-subset(data.muslim, select=c("latitude", "longitude"))
voronoi<-voronoi(df)
plot(voronoi)
proj4string(voronoi)<- CRS("+proj=longlat +datum=WGS84")

W_nb <- poly2nb(voronoi, queen = T)
##visualize queen neighbors
plot(voronoi, col="gray", border="blue")
plot(W_nb, df, add=TRUE, pch=20)
#Compute spatial weights:  A spatial weights matrix reflects the intensity of the geographic relationship between observations
W_cont_matW <- nb2listw(W_nb, style = "W", zero.policy = T)
summary(W_nb)
plot(W_nb,coordinates(df))

m2<-glm.nb(formula=number~ elite_house+ No_household, data=data.muslim)
mirr2<-negbinirr(formula=number~ elite_house+ No_household, data=data.muslim)
summary(m2)
resm2 <- m2$residuals
moran.mc(resm2, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I =0.04 indicates no spatial autocorrelation and no need for spatial filtering

m3<-glm.nb(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.muslim)
mirr3<-negbinirr(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.muslim)
summary(m3)
resm3 <- m3$residuals
moran.mc(resm3, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.03 indicates no spatial autocorrelation and no need for spatial filtering


m4<-glm.nb(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=data.muslim)
mirr4<-negbinirr(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=data.muslim)
summary(m4)
resm4 <- m4$residuals
moran.mc(resm4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.03 indicates no spatial autocorrelation and no need for spatial filtering


texreg(list(mirr2,mirr3,mirr4), stars=c(0.01, 0.05,0.1), digits=3, file="TableA10.tex")



### Figure A11###

pdf("waternetwork.pdf",   width=7, height=5)
ggmap(basemap_transparent) + geom_point(
  aes(x = longitude, y = latitude),
  data = data.nowater, colour = "red", position="jitter", alpha=1, size=2, shape=17)+geom_point(
    aes(x = longitude, y = latitude),
    data = data.water, colour = "blue", position="jitter", alpha=0.8, size=1)+
  geom_point(
    aes(x = longitude, y = latitude),
    data = palace, colour = "purple", alpha=0.8,shape=18, size=5)
dev.off()

#### Table A11 ####
df<-subset(data.water, select=c("latitude", "longitude"))
voronoi<-voronoi(df)
plot(voronoi)
proj4string(voronoi)<- CRS("+proj=longlat +datum=WGS84")

W_nb <- poly2nb(voronoi, queen = T)
##visualize queen neighbors
plot(voronoi, col="gray", border="blue")
plot(W_nb, df, add=TRUE, pch=20)
#Compute spatial weights:  A spatial weights matrix reflects the intensity of the geographic relationship between observations
W_cont_matW <- nb2listw(W_nb, style = "W", zero.policy = T)
summary(W_nb)
plot(W_nb,coordinates(df))


m1<-glm.nb(number~ NonMuslim + No_household, data=data.water)
mirr1<-negbinirr(number~ NonMuslim + No_household, data=data.water)
summary(m1)
resm1 <- m1$residuals
moran.mc(resm1, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.04 indicates no spatial autocorrelation and no need for spatial filtering

m2<-glm.nb(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
mirr2<-negbinirr(number~ NonMuslim + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
summary(m2)
resm2 <- m2$residuals
moran.mc(resm2, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.04 indicates no spatial autocorrelation and no need for spatial filtering


m4<-glm.nb(formula=number~ elite_house+ No_household, data=data.water)
mirr4<-negbinirr(formula=number~ elite_house+ No_household, data=data.water)
summary(m4)
resm4 <- m4$residuals
moran.mc(resm4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I =0.04 indicates no spatial autocorrelation and no need for spatial filtering

m5<-glm.nb(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
mirr5<-negbinirr(number~ elite_house + No_household+ Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
summary(m5)
resm5 <- m5$residuals
moran.mc(resm5, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.05 indicates no spatial autocorrelation and no need for spatial filtering


m7<-glm.nb(number~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
mirr7<-negbinirr(number~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+ Eyup, data=data.water)
resm7 <- m7$residuals
moran.mc(resm7, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I = 0.02 indicates no spatial autocorrelation and no need for spatial filtering

texreg(list(mirr1,mirr2,mirr4,mirr5,mirr7), stars=c(0.01, 0.05,0.1), digits=3,reorder.coef = c(1, 3, 2), file="tableA11.tex",omit.coef="(Elevation)|(WalkingDistanceKilometers)|(StraightLineDistanceKilomete)|(Eyup)",custom.coef.names = c("Non-Muslim Neighborhood", "No. of households","Elite neighborhood"),out="TableA11.tex")

#### Table A12 ####
ols1<-lm(water_per_day~ NonMuslim + No_household, data=data.water)
summary(ols1)
resols1 <- ols1$residuals
moran.test(resols1, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.18 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols1<-SpatialFiltering(water_per_day ~ NonMuslim + No_household, data=data.water, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
ols1Eig<-lm(water_per_day~ NonMuslim + No_household+fitted(Eigenols1), data=data.water)
summary(ols1Eig)


ols2<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
summary(ols2)
resols2 <- ols2$residuals
moran.test(resols2, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.16  indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols2<-SpatialFiltering(water_per_day ~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new 
ols2Eig<-lm(water_per_day~ NonMuslim + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(Eigenols2), data=data.water)
summary(ols2Eig)



ols4<-lm(water_per_day~ elite_house + No_household, data=data.water)
summary(ols4)
resols4 <- ols4$residuals
moran.test(resols4, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.21 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols4<-SpatialFiltering(water_per_day ~ elite_house + No_household, data=data.water, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols4Eig<-lm(water_per_day~ elite_house + No_household+fitted(Eigenols4), data=data.water)
summary(ols4Eig)


ols5<-lm(water_per_day~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
summary(ols5)
resols5 <- ols5$residuals
moran.test(resols5, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.19 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols5<-SpatialFiltering(water_per_day ~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols5Eig<-lm(water_per_day~ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+fitted(Eigenols5), data=data.water)
summary(ols5Eig)


ols7<-lm(water_per_day~ NonMuslim + elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water)
summary(ols7)
resols7 <- ols7$residuals
moran.test(resols7, listw=W_cont_matW, zero.policy=T)
#Moran I statistic of 0.17 with p=7.161e-07 indicates mild spatial autocorrelation
old <- Sys.time()
Eigenols7<-SpatialFiltering(water_per_day ~ NonMuslim+elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup, data=data.water, nb = W_nb, style="W", tol=0.2, ExactEV = T, zero.policy = T)
new <- Sys.time()    
old-new #2.96 min.
ols7Eig<-lm(water_per_day~ NonMuslim+ elite_house + No_household+Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete+Eyup+fitted(Eigenols7), data=data.water)
summary(ols7Eig)

stargazer(ols1Eig,ols2Eig,ols4Eig,ols5Eig,ols7Eig,digits=2,omit=c("Elevation","WalkingDistanceKilometers","StraightLineDistanceKilomete","fitted"), covariate.labels = c("Non-Muslim neighborhood","Elite neighborhood", "No. of households","Constant"), omit.stat=c("rsq","ser","f"), add.lines=list(c("Geographic Controls", "No", "Yes","No", "Yes", "Yes"), c("Moran Eigenvectors", "Yes", "Yes", "Yes", "Yes", "Yes")), column.labels=c("Model 1", "Model 2", "Model 3", "Model 4","Model 5"), star.cutoffs=c(0.1, 0.05, 0.01), star.char=c("*", "**", "***"), notes = c("* p<0.1; ** p<0.05; *** p<0.01"), notes.append = F,out="TableA12.tex")


#### Table A13 ####
ch1<-hurdle(number~ NonMuslim  + No_household, data=data.water,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
summary(ch1)
resch1 <- ch1$residuals
moran.mc(resch1, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.05 indicates no spatial autocorrelation

ch2<-hurdle(number~ NonMuslim  + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch2<-residuals(ch2)
moran.mc(resch2, listw=W_cont_matW, zero.policy=T,nsim=10000)
summary(ch2)
#Moran I statistic of 0.04 indicates no spatial autocorrelation

ch4<-hurdle(number~ elite_house + No_household, data=data.water,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
summary(ch4)
resch4 <- ch4$residuals
moran.mc(resch4, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.04 indicates no spatial autocorrelation


ch5<-hurdle(number~ elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch5<-residuals(ch2)
moran.mc(resch5, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.04 indicates no spatial autocorrelation
summary(ch5)


ch7<-hurdle(number~ NonMuslim+elite_house + No_household + Elevation + WalkingDistanceKilometers+ StraightLineDistanceKilomete, data=data.water,
            dist = "negbin",
            zero.dist = "binomial",
            link = "probit")
resch7<-residuals(ch7)
moran.mc(resch7, listw=W_cont_matW, zero.policy=T,nsim=10000)
#Moran I statistic of 0.02 indicates no spatial autocorrelation
summary(ch7)

texreg(list(ch1,ch2,ch4,ch5,ch7), stars=c(0.01, 0.05,0.1), digits=3,file="TableA13.tex")

